#! /usr/bin/python
# VERSION 1.0 Peter Daldrop 28/1/2008
#changelog:1.0 all basic features
import os,sys,string,math,random
import lapack_lite
from Numeric import *
from numpy import *
#from openeye.oechem import *
####Defining functions
def sphconv(file_in):

#output file
    file_out = file_in+".pdb"
    f_in = open(file_in,'r')
    f_out = open(file_out , 'w')

    start = 'false'

    atom_type_dict ={7:'C', 3:'O', 4:'N', 6:'I'}
    

    for i in f_in.xreadlines():
        if start == 'true':
            #print i
            number, x, y, z, junk1, junk2, junk3, colour = string.split(string.strip(i))
            out_line = 'ATOM' + string.rjust(number, 7) + string.rjust(atom_type_dict[int(colour)], 3) + string.rjust('SPH', 6)+ string.rjust(number,6) + string.rjust(x[0:5],12) + string.rjust(y[0:6],8) + string.rjust(z[0:6],8) 
            f_out.write(out_line + '\n')
            f_out.write('TER\n') 
      
        elif i[0:7] == 'cluster':
            start = 'true'
     
    f_out.close()
    f_in.close() 
def coordout(coord):
	out="  "+str(coord)+"00000"	
	out=string.split(out,".") #split at .
	out=out[0][-3:]+"."+out[1][:5] #refuse with right format
	return out #return correct string
def dist(x,y): #calculates distance between x and y coordinates (x and y should be lists of three coordinates)
	dist=dot((x-y),transpose((x-y)))
	return math.sqrt(dot((x-y),transpose((x-y))))

def errorinput(typ):
	print "Syntax error!"
        if typ==1:
		print "Too few arguments given!"
	if typ==2:
		print "Too many argumets given!"
	print "Syntax: boxnibbler <input spheres><receptor file><output spheres>{critical distance}"
	print "input spheres: sph file of the full sphere set"
	print "receptor file: pdb file of the docking prepared receptor containing the atoms the spheres shouldn't interfere with"
	print "output spheres: file name for the generated sphere set"
	print "critical distance: if sphere in input set below this distance from any atom of the receptor file it will be discarded (default is 2)"

### This script removes all spheres from a given sphere set clashing with protein atoms in a specified pdbfile. The resulting sphere set is then written into an output file. The openeye python library is necessary for this script.
clashdist=2 #critical distance below sphere is regarded clashing with protein atom
if sys.argv.__len__()<4:
	errorinput(1)
	sys.exit()
elif sys.argv.__len__()==4:
	print "critical distance set to default"
elif sys.argv.__len__()==5:
	clashdist=float(sys.argv[4])
	print "critical distance set to "+sys.argv[4]
elif sys.argv.__len__()>5:
	errorinput(2)
	sys.exit()
else:
	print "weired error contact developer for help"
sph_in=open(sys.argv[1],'r') #containing input sphere list
prot_in=open(sys.argv[2],'r') #containing pdb with protein atoms
out_file=open(sys.argv[3],'w') #file to put the resulting spheres
sphcoord=[] #array for the sphere coordinates
protcoord=[] #array for the protein coordinates
proxcoord=[] #array for the relevant proximal protein coordinates
validsph=[]
for i in sph_in.readlines(): #for all lines
	if i[0]==" ": #true if line contains sphere coordinates
		sphcoord.append(array((float(i[6:14]),float(i[16:24]),float(i[26:34]) ))) #adds coordinates to list with all  sphere's coordinates
	elif i[0:7]=='cluster':
		print "\nHeader succesfully read in\n" #log output	
	else: 
		out_file.write(i) #write out header line
sphcoord=array(sphcoord)
print "Sphere set was read in successfully!\n"
del i #deleting previous counter variable
for i in prot_in.readlines(): #read in all lines in pdb file
	if i[0:4]=="ATOM": #if line contains protein ATOM
		protcoord.append(array([float(i[31:38]),float(i[39:46]) ,float(i[47:54])])) #extract coordinates for atom and store in list			
protcoord=array(protcoord) #make array from list
print "Protein atoms read in successfully\n"
rim=array([4,-4]) #defines amount of space around ligand box to be used for clash analysis
#print sphcoord
xbounds=array([max(sphcoord[:, 0]),min(sphcoord[:, 0])])+rim # x coordinates range
ybounds=array([max(sphcoord[:, 1]),min(sphcoord[:, 1])])+rim # y coordinates range
zbounds=array([max(sphcoord[:, 2]),min(sphcoord[:, 2])])+rim # z coordinates range
n=0
while n<protcoord.shape[0]-1: # all protein atom coordinates 
	if (protcoord[(n),0]<xbounds[0] and protcoord[n,0]>xbounds[1]): #if in x range
		if (protcoord[(n),1]<ybounds[0] and protcoord[n,1]>ybounds[1]): #if in y range
			if (protcoord[(n),2]<zbounds[0] and protcoord[n,2]>zbounds[1]): #if in z range
				proxcoord.append(array(protcoord[n,:])) #append coordinates to list of relevant coordinates
	n=n+1
	
proxcoord=array(proxcoord)  #reshaping list to array
print str(proxcoord.shape[0]) +" Atoms found in ligand environment\n" #just a log output
n=0 #initialise sphere counter variable
while n<sphcoord.shape[0]:
	m=0 #initialise atom counter variable
	clash=False #initialise boolean if clash exists
	while m<proxcoord.shape[0]:#for all protein atoms sourrounding sphere area
		if (dist(proxcoord[m],sphcoord[n])<clashdist):  #if distance lower
			clash=True
		if clash: #if clash with one Atom occurs stop checking
			break
		
		m=m+1 #increment counter
	if (clash==False): #if sphere doesnt clash
		validsph.append(sphcoord[n]) #store validated sphere's coordinates
	n=n+1 #increment counte
validsph=array(validsph) #reformat to array
n=0
clustersize="  "+str(validsph.shape[0]) #number of spheres in set
out_file.write("cluster     1   number of spheres in cluster   "+clustersize[-3:]+"\n") #write correct cluster size line
while n<validsph.shape[0]: #for all nonclashing spheres
	number="  "+str((n+1))	#generate sphere number containing string
	out_file.write("  "+number[-3:]+" "+coordout(validsph[n,0])+" "+coordout(validsph[n,1])+" "+coordout(validsph[n,2])+"   0.500    0 0  7\n")#write sphere coordintates embedded in the other information
	n=n+1
print str(validsph.shape[0])+" spheres not clashing" #log output

out_file.close() 
sphconv(sys.argv[3])
sys.exit()




